function out = getExpectedShortRates(setupStep3,outputStep3,theta2Step3)

% Computing expected 3-months rate 3-months and 1-year ahead
model       = outputStep3.model;
posMat      = (model.matSelect == 3);
fcHorizon   = 12; 
if setupStep3.modelNum == 1
    numSim      = 1D4;  
elseif setupStep3.modelNum == 2    
    numSim      = 1D2;  
end  

% For the P dynamics we use estimates from Step 3. Note that we
% here apply the sigma from step 3
nx    = size(model.gx,2);
[h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz]   = unfoldTheta2_threshold(theta2Step3,nx);

% The shocks for forecasting the states
rng(1,'twister')
shocks        = randn(nx,numSim,fcHorizon);
shocks        = shocks-repmat(mean(shocks,2),1,numSim,1);    %ensure conditional zero mean
shocks        = shocks./repmat(std(shocks,1,2),1,numSim,1);  %ensure conditional unit variance
shocksz       = randn(numSim,fcHorizon);
shocksz       = shocksz - repmat(mean(shocksz,1),numSim,1);  %ensure conditional zero mean
shocksz       = shocksz./repmat(std(shocksz,1,1),numSim,1);  %ensure conditional unit variance
T             = size(setupStep3.econAct,1);
forecasts     = zeros(fcHorizon,T);
vectorOn      = 1;
for t=1:T
    x           = outputStep3.factors(:,t);
    z           = setupStep3.econAct(t,1);
    if vectorOn == 1
        factorSim        = nan(nx,numSim,fcHorizon+1);
        factorSim(:,:,1) = repmat(x,1,numSim);
        zSim             = nan(numSim,fcHorizon+1);
        zSim(:,1)        = z;
        for i = 1:fcHorizon
            boom  = (zSim(:,i) >= 0)';
            reces = (zSim(:,i) < 0)';
            factorSim(:,:,i+1) = repmat(boom,nx,1).*repmat(h0_1,1,numSim) + repmat(reces,nx,1).*repmat(h0_2,1,numSim) ...
                + repmat(boom,nx,1).*(hx_1*factorSim(:,:,i)) + repmat(reces,nx,1).*(hx_2*factorSim(:,:,i))...
                + sigma*shocks(:,:,i); %randn(nx,1);
            if isfield(setupStep3,'numSpanned')
                % for the spanned model, the macro variable appears last in factors
                zSim(:,i+1) = factorSim(nx,:,i);
            else
                zSim(:,i+1) = repmat(gama0,numSim,1) + gamaz*zSim(:,i) + (gamax'*factorSim(:,:,i))' + sigmaz*shocksz(:,i); %randn(1,1);
            end
            % Conditional forecast of 3 month rate
            if setupStep3.modelNum == 1
                forecasts(i,t)            = model.g0(posMat) + model.gx(posMat,:)*mean(factorSim(:,:,i+1),2);
            elseif setupStep3.modelNum == 2
                tmp =  zeros(length(model.matSelect),numSim);
                for s=1:numSim
                   tmp(:,s) = modelFit(model,factorSim(:,s,i+1));
                end
                forecasts(i,t) = mean(tmp(posMat,:),2);
            end
        end
    else
        xSim        = zeros(nx,fcHorizon+1,numSim);
        xSim(:,1,:) = repmat(x,1,numSim);
        zSim        = zeros(fcHorizon+1,numSim);
        zSim(1,:)   = repmat(z,1,numSim);
        for j = 1:fcHorizon
            for s=1:numSim
                if zSim(j,s) >= 0
                    xSim(:,j+1,s) = h0_1 + hx_1*xSim(:,j,s) + sigma*shocks(:,s,j);
                elseif zSim(j,s) < 0
                    xSim(:,j+1,s) = h0_2 + hx_2*xSim(:,j,s) + sigma*shocks(:,s,j);
                end
                if isfield(setupStep3,'numSpanned')
                    % for the spanned model, the macro variable appears last in factors
                    zSim(j+1,s) = xSim(nx,j+1,s);
                else
                    zSim(j+1,s) = gama0 + gamaz*zSim(j,s) + gamax'*xSim(:,j,s)+ sigmaz*shocksz(s,j);
                end
            end
            
            % Conditional forecast of 3 month rate
            if setupStep3.modelNum == 1
                forecasts(j,t)            = model.g0(posMat) + model.gx(posMat,:)*mean(xSim(:,1+j,:),3);
            elseif setupStep3.modelNum == 2
                tmp =  zeros(length(model.matSelect),numSim);
                for s=1:numSim
                    tmp(:,s) = modelFit(model,xSim(:,j+1,s));
                end
                forecasts(j,t) = mean(tmp(posMat,:),2);
            end
        end
    end
end
% Forecast errors
forecastErrorsOfSurvey3md = setupStep3.survey(:,2)/100- forecasts(3,:)';
forecastErrorsOfSurvey1yr = setupStep3.survey(:,3)/100- forecasts(12,:)';
forecastErrorsOfYield3md  = setupStep3.data(4:end,posMat) - forecasts(3,1:end-3)';
forecastErrorsOfYield1yr  = setupStep3.data(13:end,posMat) - forecasts(12,1:end-12)';
surveyErrors3md           = setupStep3.data(4:end,posMat) - setupStep3.survey(1:end-3,2)/100;
surveyErrors1yr           = setupStep3.data(13:end,posMat)- setupStep3.survey(1:end-12,2)/100;
% Unconditional RMSEs
RMSEofSurvey3md     = sqrt(mean(forecastErrorsOfSurvey3md(~isnan(forecastErrorsOfSurvey3md)).^2));
RMSEofSurvey1yr     = sqrt(mean(forecastErrorsOfSurvey1yr(~isnan(forecastErrorsOfSurvey3md)).^2));
RMSEofYield3md      = sqrt(mean(forecastErrorsOfYield3md.^2));
RMSEofYield1yr      = sqrt(mean(forecastErrorsOfYield1yr.^2));
RMSEsurveyErrors3md = sqrt(mean(surveyErrors3md(~isnan(surveyErrors3md)).^2));
RMSEsurveyErrors1yr = sqrt(mean(surveyErrors1yr(~isnan(surveyErrors1yr)).^2));
% Conditional RMSEs based on NBER recessions
selectRec          = setupStep3.NBERrec    .*~isnan(forecastErrorsOfSurvey3md);
selectExp          = (1-setupStep3.NBERrec).*~isnan(forecastErrorsOfSurvey3md);
RMSEofSurvey3mdRec = sqrt(mean(forecastErrorsOfSurvey3md(selectRec==1).^2));
RMSEofSurvey3mdExp = sqrt(mean(forecastErrorsOfSurvey3md(selectExp==1).^2));
RMSEofSurvey1yrRec = sqrt(mean(forecastErrorsOfSurvey1yr(selectRec==1).^2));
RMSEofSurvey1yrExp = sqrt(mean(forecastErrorsOfSurvey1yr(selectExp==1).^2));

selectRec          = setupStep3.NBERrec(4:end)    .*~isnan(forecastErrorsOfYield3md);
selectExp          = (1-setupStep3.NBERrec(4:end)).*~isnan(forecastErrorsOfYield3md);
RMSEofYield3mdRec  = sqrt(mean(forecastErrorsOfYield3md(selectRec==1).^2));
RMSEofYield3mdExp  = sqrt(mean(forecastErrorsOfYield3md(selectExp==1).^2));

selectRec          = setupStep3.NBERrec(13:end)    .*~isnan(forecastErrorsOfYield1yr);
selectExp          = (1-setupStep3.NBERrec(13:end)).*~isnan(forecastErrorsOfYield1yr);
RMSEofYield1yrRec  = sqrt(mean(forecastErrorsOfYield1yr(selectRec==1).^2));
RMSEofYield1yrExp  = sqrt(mean(forecastErrorsOfYield1yr(selectExp==1).^2));

selectRec          = setupStep3.NBERrec(4:end)    .*~isnan(surveyErrors3md);
selectExp          = (1-setupStep3.NBERrec(4:end)).*~isnan(surveyErrors3md);
RMSEsurveyErrors3mdRec = sqrt(mean(surveyErrors3md(selectRec==1).^2));
RMSEsurveyErrors3mdExp = sqrt(mean(surveyErrors3md(selectExp==1).^2));

selectRec          = setupStep3.NBERrec(13:end)    .*~isnan(surveyErrors1yr);
selectExp          = (1-setupStep3.NBERrec(13:end)).*~isnan(surveyErrors1yr);
RMSEsurveyErrors1yrRec = sqrt(mean(surveyErrors1yr(selectRec==1).^2));
RMSEsurveyErrors1yrExp = sqrt(mean(surveyErrors1yr(selectExp==1).^2));

% A figure
dates = datenum(setupStep3.timeIndex);
timeIndex = str2num(datestr(dates,10))+str2num(datestr(dates,5))/12;
figure('Name','Surveys','NumberTitle','off');
subplot(2,1,1)
hold on
plot(timeIndex,setupStep3.survey(:,2)/100,'-*')
plot(timeIndex,forecasts(3,:)','-k')
axis tight
hold off
title(['Surveys for 3-monthts rate: 3 months ahead. RMSE = ',num2str(RMSEofSurvey3md)]);
subplot(2,1,2)
hold on
plot(timeIndex,setupStep3.survey(:,3)/100,'-*')
plot(timeIndex,forecasts(12,:)','-k')
axis tight
hold off
title(['Surveys for 3-monthts rate: 1 year ahead. RMSE = ',num2str(RMSEofSurvey1yr)]);

% Output
out.forecasts                 = forecasts;
out.forecastErrorsOfSurvey3md = forecastErrorsOfSurvey3md;
out.forecastErrorsOfSurvey1yr = forecastErrorsOfSurvey1yr;
out.forecastErrorsOfYield3md  = forecastErrorsOfYield3md;
out.forecastErrorsOfYield1yr  = forecastErrorsOfYield1yr;
out.surveyErrors3md           = surveyErrors3md;
out.surveyErrors1yr           = surveyErrors1yr;
out.RMSEofSurvey3md           = RMSEofSurvey3md;
out.RMSEofSurvey1yr           = RMSEofSurvey1yr;
out.RMSEofYield3md            = RMSEofYield3md;
out.RMSEofYield1yr            = RMSEofYield1yr;
out.RMSEsurveyErrors3md       = RMSEsurveyErrors3md;
out.RMSEsurveyErrors1yr       = RMSEsurveyErrors1yr;
out.decom.RMSEofSurvey3mdRec  = RMSEofSurvey3mdRec;
out.decom.RMSEofSurvey3mdExp  = RMSEofSurvey3mdExp;
out.decom.RMSEofSurvey1yrRec  = RMSEofSurvey1yrRec;
out.decom.RMSEofSurvey1yrExp  = RMSEofSurvey1yrExp;
out.decom.RMSEofYield3mdRec   = RMSEofYield3mdRec;
out.decom.RMSEofYield3mdExp   = RMSEofYield3mdExp;
out.decom.RMSEofYield1yrRec   = RMSEofYield1yrRec;
out.decom.RMSEofYield1yrExp   = RMSEofYield1yrExp;
out.decom.RMSEsurveyErrors3mdRec = RMSEsurveyErrors3mdRec;
out.decom.RMSEsurveyErrors3mdExp = RMSEsurveyErrors3mdExp;
out.decom.RMSEsurveyErrors1yrRec = RMSEsurveyErrors1yrRec;
out.decom.RMSEsurveyErrors1yrExp = RMSEsurveyErrors1yrExp;